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Abstract. We present a complete analysis of the dynamics of a Bosc-Einstcin 
condensate trapped in a symmetric triple-well potential. Our classical analogue 
treatment, based on a time-dependent variational method using SU(3) coherent states, 
includes the parameter dependence analysis of the equilibrium points and their local 
stability, which is closely related to the condensate collective behaviour. We also 
consider the effects of off-site interactions, and how these 'cross-collisions' may become 
relevant for a large number of trapped bosons. Besides, we have shown analytically, by 
means of a simple basis transformation in the single-particle space, that an integrable 
sub-regime, known as twin-condensate dynamics, corresponds in the classical phase 
space to invariant surfaces isomorphic to the unit sphere. However, the quantum 
dynamics preserves the twin-condensate defining characteristics only partially, thus 
breaking the invariance of the associated quantum subspace. Moreover, the periodic 
geometry of the trapping potential allowed us to investigate the dynamics of finite 
angular momentum collective excitations, which can be suppressed by the emergence 
of chaos. Finally, using the generalized purity associated to the su(3) algebra, we were 
able to quantify the dynamical classicality of a quantum evolved system, as compared 
to the corresponding classical trajectory. 



PACS numbers: 03.65. Sq, 03.75.Lm, 03.75.Kk 
1. Introduction 



The experimental observation of the Bose-Einstein condensation in systems of ultracold 
dilute alkali atomic clouds pQ has led to an active development of various studies 
concerning collective quantum phenomena. The special property of coherence presented 
in condensates of a large number of atoms allowed mean-field and semiclassical 
approaches to discover the fascinating dynamics of such a large quantum system. 

The prospect of using magneto-optical traps with diverse confining geometries 
allowed the study of various quantum dynamical regimes, also including the possibility 
of testing the frontier between the quantum and classical mechanics through the increase 
of the total particle number to the macroscopic scale. 

Several surveys have been performed in the case of a Bose-Einstein condensate 
trapped in a double well-potential under the two-mode approximation. Various subjects 
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were considered in these studies, such as tunneling dynamics [2J, entanglement of modes 
[3] and quantum phase transition [I]. More recently the effects of the so-called cross- 
collisions have also been considered, and it was shown that these usually neglected 
higher order approximation terms can become quite important for large numbers of 
trapped bosons [5]. 

The dynamics of condensates in a triple- well potential has already been investigated 
a few times in the literature [6]. However, little attention has been paid to the specific 
case in which all three wells are identically coupled in pairs [7]. This arrangement is very 
interesting because it represents the simplest trapping potential where we can observe 
the rotation of the condensate. 

In the present work we examine the dynamics of a condensate in a symmetric triple- 
well potential, under the three-mode approximation and considering the cross-collision 
terms. First, we investigate the classical analogue of the model and its phase space 
structure, by determining its fixed points and the Hamiltonian flow in the vicinity of 
these special solutions, which is closely related to the condensate collective excitations. 
Then, we choose the coherent states located in the proximity of the classical stationary 
states as initial quantum states. This enables us to capture in the quantum dynamics 
the characteristics associated to each of the classical fixed points. This strategy has 
been used before to observe the effect of chaos in the entanglement dynamics [8]. 

In section [2J we derive the model in the three-mode approximation, including the 
cross-collisional terms, and the corresponding Hamiltonian in terms of the SU(3) group 
generators. We also present the classical approximation based on the time- dependent 
variational principle (TDVP) [9] and the SU(3) coherent states [10]. The equilibrium 
points of the classical analogue equations of motion are divided into three classes: twin- 
condensate states, single depleted well states and vortex states. Here we show in detail 
the stability diagrams of the fixed points as functions of two collision parameters. The 
behaviour of the lowest energy fixed point and the occurrence of a quantum phase 
transition has been reported recently [TTj . 

In section [3] we relate a simple permutation symmetry of the Hamiltonian to a 
transformation of the single-particle basis that leads to the twin-condensate subspace. 
Then, we demonstrate that the classically invariant subspace of twin-condensates is not 
fully preserved quantum mechanically. In addition, we compare the quantum population 
dynamics in the local modes with the corresponding classical evolution. Particularly 
interesting is the behaviour of the vortex states, for which we show the dynamics of a 
quantity associated with the condensate angular momentum along the symmetry axis 
of the trapping potential. The variation of the collision parameters can lead the vortex 
solutions to chaotic dynamics, and consequently inhibiting the condensate rotation in a 
preferred sense. It is notable that this model present chaotic behaviour at the classical 
level. Moreover, similarly to previous results for the double-well condensate model p[| 
with the measure known as generalized purity [12] , we make a quantitative comparison 
between the quantum evolution and the classical approximation based on the coherent 
states. Section [4] is reserved to our conclusions and final considerations. 
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2. Model 

Suppose a spinless boson of mass m trapped in a potential V(f), whose equivalent 
minima are located at fj = (xj, yj, 0), j = 1, 2, 3: 

V W = ^ fl [(* - + (v- mf] + (2-i) 

Note that the potential is harmonic in ^-direction with angular frequency u and 
we choose V(fj) = 0. The parameter q represents the distance between the potential 
minima and the origin of the coordinate system. For simplicity, we also choose V(r) 

(2) —* 

so that the harmonic approximation V- if) around each minimum is isotropic with 
angular frequency to. We can assign to each potential (r), taken independently, a 
localized ground state: 

(x - Xjf + {y- yj) 2 + z 2 



(f]uj) =Uj{r) = ^-yexp 



2d 2 



(2.2) 



The Gaussian width is described by the parameter d = ^Jh/muj. Considering that 
q 3> d, the states \uj) are almost orthogonal, since e = (uj\uk) = exp [— 3q 2 /4d 2 } <C 1, 
for j ^ k. 

The states \uj) are not eigenstates of the full single-particle Hamiltonian H = 
p 2 /2m + V(r), but these localized states generate the same eigenspace of the three 
lowest energy states of H at first order in e. That is, we have the following single- 
particle lowest energy states at 0(e): 

r 1 

|e 2 ) = -^(K) + K)-2K)) (2.3) 

|e 3 ) = - l n 2)) 

The states \e<i) and \e-&) are degenerated and the energy of these states differs from 
the energy of \e\) by 3f2, with the tunneling rate defined by: 

= ( Uj \V(f) - V} 2 \f)\u k )+e^- j + k. (2.4) 

Generally, for q ^> d, fl is negative and \ei) is the single-particle ground state, 
whose energy difference for the next two states is proportional to e. 

So far we studied only the single-particle potential. Now, considering a system of 
iV condensate bosons, we introduce the many-particle Hamiltonian: 

H = J d 3 r^{r)Hij;(r) + ^ J d 3 r d 3 r'ft(r)ft(r' )V(f, r )^(f). (2.5) 

where H is the single-particle Hamiltonian and V(r, r 1 ) is the pair interaction potential. 
Only pair interactions are significant, because we assume dilute condensate gas, 
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characterized by a high mean free path, which makes higher order collisions unlikely. 
For particles of momentum approaching zero, as in the case of a condensate, a good 
approximation for the interaction potential is the following effective potential [13] : 

V(f, 7*) = Z^8(f- r) = V S(f- r). (2.6) 
m 

In the previous equation we defined the s-wave scattering length, denoted by a. 
Now, supposing that the states of (12.31) are the only significantly populated, we get our 
second approximation in (I2.5p . namely, the three-mode approximation: 

-0(r) ps (r\u 1 )a 1 + (f\u 2 )a 2 + (f\u 3 )a 3 ; (2.7) 

where we defined a,j as the annihilation operator for the state \uj). Therefore, the field 
operator ip{r) may be expanded in the local modes described in ( I2.2p . which generate 
the same single-particle space as the energy states of ( 12. 3p . Substituting ( 12. 6 j) and ( 12. 7p 
in (I2.5p . we obtain the following Hamiltoniarjjl, retaining only the terms up to order e§: 

H = Q' ajafc + k a j 2fl j — 2A ajaja|.a m . (2.8) 

j^k j j,k,m 

The third summation symbol in ( 12. 8 j) indicates the sum over three different indices. 
In the previous equation we defined the effective tunneling rate f2' = Q + 2A(A r — 1) 
and the following collision parameters: 



^0 /"^3 4 /-a 
k = — / d r uAr) 



A = — / d 3 r u^{r)uk{r) = ke" 2 ] j ^ k. 



(2.9) 



2 i J 227T2<i 3 

The self- collision rate k is proportional to the collision frequency between bosons in 
the same potential well, while the cross-collision rate A is proportional to the collision 
frequency between bosons coming from different local modes. Note that the effective 
tunneling rate Q! depends on the cross-collision parameter and the total number of 
trapped bosons. Also, observe that we have |A| <C |ac| for e <§C 1. Although the cross- 
collision rate has a lower order of magnitude than Q, the product A^A in Q' may be 
relevant for large numbers of trapped bosons. Therefore, the cross-collisions become 
important in comparing theoretical and experimental results, because they include 
significant effects of the number of particles in the effective tunneling rate. 

Since H conserves the total number of trapped particles under unitary evolution, 
we can use the following homomorphism between the bilinear bosonic operators and the 
generators of the SU(3) group: 

Qi = ^i a \ a i - 4^), Q2 = g(4°i + 4 a 2 - 24*33), ^ 
j k = i(4«i - 4 Gfc )' p k = a l a j + a ) a k] 

\ Note that we use the bosonic canonical commutation relations to simplify the final expression for H. 
Wc also discard the terms dependent only on the conserved total number of particles, as such terms do 
not change the system dynamics for a fixed value of TV. 
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for k = 1,2,3 and j = (k + l)mod3 + 1. The Hamiltonian f!2.8p can be rewritten 
in terms of the operators (I2.10p . characterizing SU(3) as the dynamical group of the 
systeirj§|: 



H=\tt- 2Ajj (Pi + P 2 + Ps) + '-(AQ\ + 3Ql) 

+ h[2Q x {P x - P 3 ) + Q 2 (2P 2 - P x - P 3 )}. 

Notice that the tunneling term is linear in the SU(3) generators, while the collision 
terms are quadratic. Within the three-mode approximation, a single-particle operator 
A and its corresponding many-particle operator A are related by: 

3 

A=^2( Uj \A\u k )a]a k . (2.12) 
j,k=i 

Therefore, considering terms up to first order in e, we can find relations between 
the generators of SU(3) and basic system observables. First, we have the following 
condensate position operators, obtained by choosing r x = — |g 0) ^Qo, OJ , r 2 = 

(-|?o, -^Qo, o) and f 3 = (q , 0, 0): 



3q eq 
x = -—Q 2 + — 



(2.13) 



4 

The operators J k , k = 1,2,3, generate a subalgebra of su(3) isomorphic to su(2), 
the angular momentum algebra. Nevertheless, linear combinations of these operators 
are directly proportional to the linear momentum operators of the condensate: 

Px = TkT ( J 3 " Ji) ; 

7- (2.14) 

Py = ~W ^ 1 3 ~ 2 > ' 
However, the condensate angular momentum operator along the z-axis is 
proportional to the sum Js = J\ + -h + Jz'- 

2.1. Classical Approximation 



J s . (2.15) 



Mean field theories and classical trajectories are used widely in the literature for the 
treatment of condensate dynamics. Here we introduce our classical approach to the 
problem based on the time- dependent variational principle. Considering independent 
variations in an arbitrary state of Hilbert space and its conjugate it is known 

§ Here again we discard some terms dependent only on the fixed total number of particles N. 
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that the extremization of the following action functional is equivalent to the Schrodinger 
equatioijjj]: 



5 



dt 



1 (^) - 

2 (?/#> 



\H\ 



<V#> 



(2.16) 



We can make approximations to the Schrodinger equation by a suitable 
parametrization of the time dependence of the state restricting the state evolution 
to a subspace of the Hilbert space. The classical approximation is achieved when the 
state is restricted to the nonlinear subspace constituted only by the coherent states 
associated with the dynamical group For the SU(3) bosonic representations, we 
have the following coherent states [TO] : 



|iV; 101,102} 



ni+n2+ri3=N 



ni\n 2 \n^,\ 



Wi 1 i02 2 |ni,ra2,n 3 ) 



fN\ 



W\a\ 



w 2 a\ + a\ 



v^oTpTKfTi 



(|wi| 2 + \w 2 \ 

N 

|0>. 



+ 1) 



(2.17) 



The variables Wi,w 2 G C parametrize the subspace composed of coherent states 
for a fixed number particles N. The states {\rix, n 2 , n-s)} form the basis of the three- 
mode bosonic Fock space, where rij G IN is the occupation number associated with the 
local state \uj). Applying the variational principle to the action functional (I2.16P with 
|^) = \N;w\, w 2 ), we obtain the classical equations of motion for the complex variables 



Wf 



1W4 



\wi\ 2 + |u> 2 | 2 + 1 
N 



WjW k 



&H 

: dw* k 



(2.18) 



for j, k = 1,2 and k 7^ j. The quantity %{w\,W2,Wi,w 2 ) is the effective classical 
Hamiltonian, given by the average of H in the coherent states: 



N ~ 



n 



-2/x 



2/x) 



N 

Wlw 2 + W2W1 



Wl + W\ + W 2 + Wn 



1 + kod 



\w 2 \ 



X 





l 4 + 


w 2 


I 4 + 1) 


(kil 


2 + l 


w 2 \ 


2 + 1) 2 



(2.19) 



[|wi| 2 (w2 + W 2 ) + \w 2 \ 2 (wl + Wi) + + 



(\ Wl \ 2 + \W 2 \ 2 + l) 2 

In the previous equation we introduced new collision parameters, which allowed 
us to eliminate the dependence on iV from the average energy per particle and also 
prevented us from problems of divergence in the macroscopic limit iV — > 00. The 
quantities x an d /-t as functions of the old collision parameters are given by: 

re(JV-l) A(JV-l) 



X 



n 



n 



(2.20) 



In what follows we make h = 1 . 
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Note that \x\ 3> \fJ>\, because supposedly we have |ac| ^> |A|. Substituting the 
Hamiltonian (12. 19[) into (12. 18ft . we obtain the classical equations of motion for the three- 
mode condensate: 

( w(\w-\ 2 - 1) 

iwj =QU1 + 2/*)(l - Wj ){ Wj + w k + l) + |2 ' 

(2.21) 

■ 

-2/jl 



[(1 - Wj)w%(Wj +W k + WjWk) 


- (1 + w j )w k (\w j 


2 -l)] 






2 + 




2 + l 



again for j, k — 1, 2 and k ^ j. The TDVP based classical approximation with coherent 
states is exact for any linear Hamiltonian in the dynamical group generators, i.e., 
the classical equations of motion (I2.2ip reproduce the correct quantum evolution of a 
coherent state in the case of non-interacting bosons, for x = /i = 0. The approximation 
also becomes exact in the classical limit, which coincides with the macroscopic limit 
N — y oo, since the total number of trapped bosons plays a role equivalent to the 
reciprocal of h [H]. However, considering x an d A* as independent parameters, we 
see that the equations of motion ( 12.211) and the Hamiltonian per particle (12.191) are 
independent of N, and therefore represent the exact results in the classical-macroscopic 
limit. 

The equations of motion (12. 18ft assume the canonical form under the following 
change of dynamical variables: 

"* *J ' e "'* ; (2 - 22) 

for j = 1,2. The canonical variables Kj and 4>j have an useful and immediate physical 
interpretation. The variable Kj is the mean bosonic occupation in the j-th potential 
well, and ( N — K\ — K 2 ) is the mean population in the third well due to the conservation 
of the total number of particles. The angular variable <pj represents the phase difference 
between the portions of the condensate located at the j'-th and third wells. 



2.2. Equilibrium Points 

2.2.1. Twin- condensates According to the equations of motion ( 12.211) . the classical 
phase space of the condensate has three invariant subspaces arising from the discrete 
rotational symmetry of the trapping potential, which are described by the following 
conditions: 

wi = W2, wi = 1, W2 = 1. (2.23) 

Each one of these equivalent invariant subspaces follows from the equality in phase 
and mean occupation of a pair of local modes, according to the transformation (I2.22p . 
Therefore, if two localized condensates have the same initial classical state, so they 
remain identical during their classical evolution. This integrable sub-regime of the 
model is known as twin- condensate dynamics. 
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The majority of the equilibrium points of the equations of motion (12.211) are located 
in the twin-condensate subspaces under the additional constraint u>i, u>2 £ K» Therefore, 
by making w± = W2 € R and Wj = in ( 12.21ft . we obtain the following polynomial 
equation for the location of the twin-condensate fixed points: 

4(1 + n)w\ - 2(1 + x + ±yi)w\ + 6fiwl + (2x - l)wi - 1 - 2/x = 0. (2.24) 

The above equation clearly has the solution w\ — W2 — 1, regardless of the choice of 
the collision parameter values. This equilibrium point, according to the interpretation 
of the transformation ( 12. 22ft . corresponds to the state with identical phase and mean 
occupation in all three local modes. We denote this solution by s±, which is the 
intersection of the three invariant surfaces in (I2.23p . Excluding the solution s\ from 
(12.24ft . we obtain a cubic equation for the remaining real fixed points: 

1 2(1 +fj) 1 2(1 + fx) 4(1 + //) V 7 

The number of real solutions of the equation ( I2.25P depends on the sign of the 
following discriminant: 

A = - X 4 - 2(3 + 7/i) X 3 + 3(2 - ll/i 2 ) X 2 + 2(5 + 12/x - 18/i 2 - 52/i 3 ) X 2 
+ 2(9 + 76// + 228/i 2 + 264/i 3 + 76/i 4 ). 

This discriminant has two real roots in x, the positive (negative) root is denoted 
by x+(/-0 (X-(/ i ))- W X-(/-0 < X < X+(a*)) the equation ( 12.25ft has only one real solution 
that we denote by s^ix^lA- However, for x — X±(a*) we have a new equilibrium point, 
which undergoes a bifurcation for x > X+(/ i ) an d X < X-(/-0> giving rise to the solutions 
denoted by s 3 (x, a0 and s 4 (x, A*)- Figured] shows the number of real solutions of equation 
(I2.24p as a function of the collision parameters. Note that for x = | + A 4 we have only 
three real solutions to (12.241) . because on this curve the solutions s\ and sz(x, A*) become 
degenerate. 

By the linearization of the equations of motion ( 12.21ft . we can analyze the stability 
of the dynamics in the vicinity of each equilibrium point. Although the position of 
the equilibrium point Si does not depend on the collision parameters, we can show the 
following condition of stability for this equilibrium point: 



in 



2 



K = ^(3 + 4aO(4 X - 9 - 4/i) < 0. (2.27) 

Therefore, considering ll > — |, si is stable for x < 1 + A*- Note that this critical 
curve of stability change coincides with the curve of degeneracy of Si and s 3 (x,a0- 

The behaviour of the fixed point S2(x, A 4 ) as a function of the collision parameters 
is shown in figure [2j(a). Therefore, for \x\ S> |a*|, S2 is generally stable (unstable) for 
X < (x > 0). Also for \x\ we can demonstrate that the equilibrium point S3 (34) 

is unstable (stable) in its region of existence in the parameter space, thus characterizing 
a saddle-node bifurcation^. 

% Wc do not show the algebraic conditions of stability for S2, S3 and S4 because they are exceedingly 
long and complicated, due mostly to the position dependence of these equilibrium points on x an d fx- 




Figure 1. [colour] Number of real solutions of equation (|2.24[) as a function of the 
collision parameters. The light grey [green] curves represent the transitional situation 
of existence of only three real solutions. Note that only the quadrants where \[l > 
are physically accessible, since the two collision parameters must have the same sign. 
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Figure 2. [colour] (a) Stability analysis of the equilibrium point S2(x, fi) as a function 
of collision parameters, (b) Stability analysis of the single depleted well states. Note 
that these are unstable equilibrium points on the curve \ — M- 



2.2.2. Single Depleted Well States Returning to the equations of motion, we can 
find two more real equilibrium points with coordinates independent of the collision 
parameters: 

Wj = -1, w k = 0; (2.28) 

for j, k = 1, 2 and j ^ k. The fixed points above are known as single depleted well 
states, and according to the transformation (I2.22p they always possess one completely 
empty local mode, while the two other modes have opposite phases and the same average 
populationEl The single depleted well states have the following stability conditions, also 

+ The rotational symmetry of the trapping potential indicates the presence of a third single depleted 
well state, corresponding to the complete depiction of the third well. However, according to the 
transformations (|2.22[) . this third equilibrium point would receive the coordinates wi,W2 — s- oo. 
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shown in figure [21(b): 

>?sdw,± =^r {-[9 + 4/i(10 + 11/i) + X (2 + x)} ± [(9 + 16/i - X ) x 

2 i (2.29) 

(9 + 128^ 3 - 8/i 2 ( X - 19) + X (5 + 3 X - X 2 ) + 16^(4 X + X 2 ))] 1 } < 0. 

2.2.3. Vortex States Within the classical approximation, the angular momentum of 
the condensate along the symmetry axis of the trapping potential is proportional to: 

2N 

(N; w ± , w 2 \J s \N; w u w 2 ) = --— — — Im (iwi - w 2 + wlw 2 ) ■ (2.30) 

\Wi\ z + \W2\ + 1 

Therefore, the classical angular momentum is directly related to the imaginary 
parts of the variables w\ and w 2 . Since all equilibrium points found previously have 
real coordinates, they represent irrotational condensate configurations. However, the 
equations ( 12. 2 1 f) have a further pair of fixed points: 

Wl = e M3 £, W2 = e^. (2.31) 

According to (]2.22p . these latter equilibrium points have the same mean occupation 
in all three local modes, since \wi\ = \w2\ = 1, but the phase difference between each pair 
of local condensates is which is the angle of rotational symmetry of the trapping 
potential. The states corresponding to (12.311) are known as vortex states, due to their 
nonzero angular momentum, which is proportional to the total number of condensate 
bosons: 

(N;w! =w* = (P^lJslN'tW! = w* 2 = e^) = ±V%N. (2.32) 

Note that the two vortex states are equivalent, because they differ only in their 
sense of rotation. The stability condition for the vortex states is given by: 

X l,± h 10 V " 16 M 6 + X) - 3(9 + A X ) 

(2.33) 

< 0. 



±^3(3 + 4//) (3 + 8//) 2 (9 + 4/i + 8*) 

Therefore, considering \i > — |, the vortex states are stable for x > ~{Aji + 9)/8, 
except on the curve x = 2/x(6 + ll/x)/(3 + 4/i). 



3. Condensate Dynamics 

3.1. Twin-condensate Dynamics 

The integrable sub-regime of twin-condensates is a direct consequence of the symmetry 
of H under the permutation of indices of the three local modes. The equivalence of 
the three local modes remains in the Hamiltonian ( 12. 19ft . which has the permutation 
symmetry of the quantities awi, aw 2 and a, where a = (|wi| 2 + |w 2 | 2 + l) 1 ■ 
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Without loss of generality, we study the twin-condensate sub-regime under the 
restriction w\ — w 2 , because the invariant surfaces in ( I2.23P are dynamically equivalent. 
Applying this condition to the coherent states in ( 12. 17ft . we obtain: 

.-1 N 



\N; wi = w 2 ) 



E 



Wia\ + wia\ + a\ 

y/2\ Wl \ 2 + l 

\/M (iuia{ + 



10} 



1112 



(3-1) 



m 1 !m 2 ! 



(2\wi\ 



N 
2 



■|0>. 



m\+m,2=N 

Now it is convenient to change the basis of the single-particle Hilbert space. Such 
transformation may be described by the following linear combinations of the bosonic 
operators: 



~7T 



V2 



(3.2) 



The operator b\ (b\) creates a particle in the equiprobable superposition of the local 
states \ui) and \u 2 ) arranged with identical (opposite) phases. So b\ is responsible for 
the occupation of the twin-condensates, while b\ represents the solitary mode, which is 
characterized by the single-particle state \u 3 ). Note that b\ is also the creation operator 
in the state |ea), found in (12. 3p . Applying the transformation (13. 2p in ( 13. ip . we have: 



\N;wi = w 2 ) 



E 

mi+m2=N 







\mi,m 2 ,0) 


ymi!m 2 ! J 


(2 




2 + l)f 



(y/2 Wl )b\ + b\ 
J2\ WI \ 2 + 1 



-1 N 



(3.3) 



|0) = \N;V2 Wl ) 



su(2)- 



The set of states {\mi, m 2 , m^)} constitutes a new basis for the bosonic Fock space, 
where rrij represents the eigenvalue of the number operator bjbj. Note that the SU(3) 
coherent state restricted to the subspace twin-condensates has no occupation in the mode 
related to the operator 6363, also called the opposite phase mode, evidencing the phase 
equality between the twin-condensates in the classical dynamics. Besides, the coherent 
state ( 13. 31) is also a SU(2) coherent state with complex coordinate y/2w± = e~ lfli tan|, 
indicating that the invariant surface of twin-condensates is topologically isomorphic to 
the unit sphere (S 2 ) [T5j . 

The states with zero mean occupation in the opposite phase mode, given by 



3"3| 



or 



E 



mi, m 2 , 0) with c miim2 G C, satisfy the following 



m\+rri2=N 
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identities: 



\a[ai + 1 
\a\at\ij;) = (7 



\a\a 2 + a\a\\ 



\a\a 2 \ 



(3.4) 



0. 



Therefore, the quantum dynamics preserves the population equality between the 
twin-condensates, because the average of the population imbalance operator Qi remains 
zero during the quantum evolution for all the initial states in the subspace generated by 
the basis {\mi, m 2 , 0)}. However, the quantum dynamics does not preserve the identical 
phases of twin-condensates, unlike the classical dynamics, because the average of b% 
does not remain constant, in spite of its initial zero value. 

At this point it is interesting to define the generators of another subalgebra of su(3) 
isomorphic to su(2), but related to the modes present in the twin-condensate sub- regime: 



b\b 2 + b\h _ .b\h-b\b 2 _ b\h-b\b 2 



S a 



s z 



(3.5) 



2 ' y 2 

Now, employing the identities f)3.4j) . we can show that: 

{i/j\S z \ip) = -(ip\a{ax + a\a 2 - a\a^\i>). 

Thus, S z is the population imbalance operator between the twin-condensates and 
the solitary mode, but only if the average occupation of the opposite phase mode is 
zero. 



(3.6) 



N 0- 




N 0- 




N 0- 




Figure 3. [colour] Classical trajectories in the invariant subspace of twin-condcnsatcs 



Wi = W 2 



-^c ' tan I for N = 30, |f2| = 1, fj, = and self-collision parameter values: 
(a) \ = 1-5, (b) x = 1-98 and (c) x — 3- The [yellow] circles ([red] crosses) indicate the 
location of stable (unstable) hxcd points with respect to the complete four-dimensional 
phase space of the three-mode condensate. The angles </> and 8 are similar to the usual 
spherical coordinate angles, except that, for convenience, we take the origin of 8 along 
the negative z-axis. 



Figure [3] displays the classical dynamics of the three-mode model in the invariant 
subspace of twin-condensates, which is represented by the unit sphere, for various values 
of the self-collision rate X: but neglecting the presence of cross-collisions. Observe that 
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the Cartesian coordinates directly represent the normalized classical averages of the 



SU(2) generators: 












f 

4 = 




= w 2 \S X \N; wi 


= w 2 ) = sin 9 cos </> 




< 


h = 




= W 2 \Sy\N] W\ 


= w 2 ) = sin 9 sin 


(3.7) 




h = 

< 




= w 2 \S z \N;w 1 


= w 2 ) = — cos 9 





Figure 0(a), for % = 1.5 < x+(/ i = 0) ~ 1.971, shows the classical dynamics 
in the absence of the equilibrium points S3 and S4. Note that all orbits surround the 
stable equilibrium point s\ and, therefore, the population imbalance varies around the 
value I z (wi = w 2 = 1) = |, which represents an identical average population in all 
three local modes. Thus, there is no preferential occupation of any local mode, and this 
behaviour characterizes the dynamical regime known as Josephson oscillation (JO). In 
figure [31(b) we show the classical dynamics for x — 1-98, i.e., a self-collision rate value 
slightly higher than the critical bifurcation parameter. The bifurcation is accompanied 
by the appearance of a separatrix, which crosses the unstable equilibrium point S3. The 
separatrix is the boundary between the orbits of the JO regime and the trajectories 
around the new stable equilibrium point S4. The trajectories near S4 vary around 
negative values of I z , indicating a preferential occupation of the third well. Therefore, 
the tunneling between twin-condensates and solitary mode is suppressed, giving rise to 
the effect known as macroscopic self-trapping (MST). The JO and MST sub-regimes 
are also present in the double-well model, whose classical approximation with SU(2) 
coherent states is quite similar to the dynamics of twin-condensates [T6] . 

Figure [31(c) displays the dynamics of twin-condensates for % = 3, a parameter value 
significantly higher than x+- F° r X > f the fixed point Si becomes unstable and takes 
the place of S3 in the separatrix. Although the orbits around S3 in the subspace of 
twin-condensates are regular, we should note that this equilibrium point is not a stable 
center in the complete phase space, since S3 behaves like a saddle point in the directions 
orthogonal to the surface u>i = w 2 . Therefore, the self-trapping dynamics which favors 
the occupation of the twin-condensates, characterized by the regular trajectories with 
I z > I in the vicinity of S3, is restricted to the invariant surface. Besides, we observe 
that an increase in self-collision rate is accompanied by an expansion of the region on 
S 2 occupied by MST orbits. 

Figure [4] shows the effect of cross-collisions on the dynamics of the mean population 
imbalance I z for the initial coherent state \N, w\ = w 2 = 0) = |0, 0, N), corresponding to 
the south pole in figure [3j The solid [blue] (dashed [green]) curve in figure [H (a) displays 
the classical time evolution of I z for N = 30, \Q\ = 1, x — 4 and \i = ^ (/•* = 10 )■ 
Notice that the greater value of the cross-collision parameter inhibits the MST, due to 
its correspondent increase of the effective tunneling rate Q' = Q(l + 2/i). Therefore, the 
condensate dynamics can be completely changed by the cross-collisions, even when the 
collision parameters /1 and x have different magnitudes. 
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Figure 4. [colour] Left (right) figure shows the classical (quantum) time evolution of 
the population imbalance I z for the initial state \N, w\ — W2 = 0) = |0, 0, N), The 
solid [blue] (dashed [green]) curve displays the results for N = 30, |fi| = 1, x = 4 and 

t l = Too = u))- 



Figure HJ(b) exhibits the quantum time evolution of I z for a direct comparison 
with the classical results. Note that, unlike previous classical results, the amplitude of 
quantum oscillations does not remain constant during the evolution of the condensate, 
evidencing a quantitative disagreement between the two approaches. However, the solid 
[blue] (dashed [green]) curve properly characterizes the MST (JO) regime, because the 
fluctuations in the population imbalance remain around I z < (I z = |). Consequently, 
the qualitative agreement between the classical and quantum approaches persists, even 
for a number of particles as small as N = 30. 

3.1.1. Generalized Purity associated with the su(3) algebra The classical equations of 
motion ( I2.2ip do not depend on the total number of trapped bosons, since they represent 
the dynamics in the macroscopic limit N — > oo. Therefore, the qualitative agreement 
between the two approaches would not be expected for a small number of particles such 
as iV = 30, unlike observed in figure HI Although mean field and classical theories are 
quite common in the treatment of condensate dynamics, little is known about the quality 
of these approximations with respect to the exact quantum calculation for a microscopic 
or mesoscopic condensate. The classical dynamics is almost exact for a macroscopic 
number of particles, but for a mesoscopic or microscopic number of bosons we must 
be able to evaluate quantitatively the quality of our approximations as a function of N 
and the propagation time, because in an usual experiment the number of condensate 
particles can range from a few hundred to the order of 10 10 bosons [T7] . 

Our classical approximation consists of restricting the system time evolution to 
the nonlinear subspace of coherent states. However, the exact quantum dynamics for 
few particles can promote the departure of an initial coherent state from the classical 
subspace, thus introducing quantitative errors in our approximation. The generalized 
purity associated with the su(3) algebra, which is a measure capable of quantifying the 



Bose-Einstein condensate in a triple-well trap 



15 



proximity of a given iV-particle state 




to the classical subspace, is given by 

\2 JL /„/,|p.l„/,\2 JL /„/,! 7.U/A2' 



E 



12 



E 



J; 



12 



(3.* 



The generalized purity is limited to the interval < V su (3)(\ip)) < 1, but we have 
T y su(3){\' l l J )) — 1 if and only if is a coherent state. On the other hand, the value of 
7 : su(3)(|' ? A)) is decreasing with the "distance" of \ip) to the subspace of coherent states. 
Therefore, "P su (3)(|' ? / ; )) is a classicality measure, because the purity value is increasing 
with the classical character of taking the coherent states as the most classical pure 
states, due to their minimum uncertainty on the phase space. 
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Figure 5. [colour] Left (right) figure shows the time evolution of 7 5 su (3) ((6363)) for 
an immediate comparison with the figure |4] The solid [blue] (dashed [green]) curve 
shows the results for TV = 30, = 1, x — 4, \i = j^q (fi = jg) and initial state 
\N,wi = w 2 = 0) = |0,0, AO- 



Figure [5j(a) displays the evolution of P su (3) for the same initial state and parameters 
of figure HI Note that the state initially suffers a sudden purity loss, which is responsible 
for the quantitative disagreement between the classical and quantum results. Therefore, 
only for a very short time interval the classical approximation roughly coincides with the 
quantum dynamics, when considering a small number of particles. However, after a brief 
initial period of purity loss, V su (s\ only exhibits fluctuations around a stable value. The 
purity loss in the MST regime is inferior to the JO regime, indicating the better quality 
of the classical approximation in the self-trapping dynamics. The orbits associated with 
the JO regime traverse a larger region of the phase space compared with the localized 
MST trajectories, as shown in figure El Therefore, the quantum states related to the JO 
regime have greater delocalization (uncertainty) on the phase space, which is responsible 
for their accentuated classicality loss. 

Within the classical approximation, the invariant subspaces of twin-condensates 
are defined by the equality of phase and mean occupation between two local modes. 
Although the quantum dynamics preserves the population equality between identical 
local condensates, according to the equation f !3.4j) . the phase difference does not remain 
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zero. In general, using the transformation (13. 2 p in the Hamiltonian ( 12. 8p . we can show 
that: 



Except for x = /i = 0orA r — >• oo, when the classical approximation is exact. 
Therefore, the subspaces with no phase difference between two local modes are not 
quantum invariant. The population dynamics of the opposite phase mode for the states 
\N,wi = W2) can not be described classically, and therefore it represents a process 
of purity (classicality) loss. Figure 13(b) shows the behaviour of (b\b 3 ) for the same 
initial state and parameters of the previous figure. In comparison with the figure [5j (a), 
we observe that a higher occupation of the opposite phase mode is accompanied by 
a greater purity loss. Therefore, the nonzero occupation of the opposite phase mode 
is partially responsible for breaking the quantum-classical correspondence in the sub- 
regime of twin-condensates. Moreover, the quantity (6363) can be used as a measure of 
validity (quality) of the classical approximation in the invariant subspaces. 

3.2. Single Depleted Well Dynamics 

The sub-regime of twin-condensates is very similar to the dynamics of a condensate in 
a double-well potential, since both models have several features in common, such as the 
integrability, the SU(2) coherent states, and especially the dynamical transition from JO 
to MST. Therefore, the single depleted well states (SDWS) represent the first example 
of an entirely new population dynamics for the three-mode model. 

According to the transformation (12.221) . the SDWS have a completely empty local 
mode, while the two other modes remain with the same mean occupation and opposite 
phases. Similarly to the self-trapping regime, the SDWS represent a dynamical effect 
of tunneling suppression, but the opposite phases between the preferentially occupied 
local modes distinguish the latter regime. 

Due to the rotational symmetry of the trapping potential, all the SDWS are 
equivalent. Thus, without loss of generality, we restrict our results to the region of 
phase space in the vicinity of {wi — — 1, w<i — 0), where the second local mode is empty. 

The SDWS have fixed positions in the phase space, but their stability varies as a 
complicated function of the collision parameters x an d At, according to the condition 
(I2.29P and the figure [21(b). Figure El (a) shows the Poincare section at 02 = for \ = 5 
and = Tj^, considering only the trajectories with the same energy of the SDWS. 
Note that the section presents the stable SDWS located in the center of a bundle of 
regular orbits. The energetically accessible region of phase space is filled by orbits 
with a behaviour similar to the SDWS, since they exhibit only small fluctuations in the 
average occupation of the second local mode, while the other two modes oscillate in 
opposite phase around y. 

Figure |6j(b) shows the Poincare section in the vicinity of the SDWS for x — ~5 and 
/1 = -^r. Here we see an unstable SDWS inserted in a chaotic region of phase space. We 



d_ 

7ft 



(N; W! = w 2 \blb 3 \N; Wl = w 2 ) ^ 0. 



(3.9) 
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Figure 6. [colour] Left (right) figure shows the Poincare section at <p2 — for N = 30, 
= 1, x = 5 {x = — 5), /i = j^p, considering only those trajectories with the same 
energy of the fixed point {w\ = —1,^2 = 0). The [yellow] circle ([red] cross) shows the 
position of the stable (unstable) equilibrium point in the section. 



also observe that the SDWS is the boundary between two distinct sets of chaotic orbits, 
because the trajectories above (below) the equilibrium point in this section preserve 
the condition K\ > — (K\ < y) during their time evolution. Unlike the stable case, 
the energetically accessible region of phase space in the vicinity of the fixed point does 
not show a behaviour similar to the SDWS, since there is no persistent depletion of the 
second local mode. 
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Figure 7. [colour] Left (right) figure shows the classical population dynamics for an 
initial condition near one of the SDWS, considering N = 30, |fi| = 1, % = 5 (x = —5), 



^ — 100 " 



The regular population dynamics in the vicinity of the SDWS for x — 5 and /i = 
is illustrated in figure 0(a), employing the initial condition indicated by the square 
[green] marker in figure [6j(a). As expected, the second local mode remains almost 
empty, while the two remaining modes exhibit periodic population inversions around 
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-5 and fi = given 

an initial condition very close to the SDWS with K\ > y. Note that the second local 
mode does not remain empty, and it presents successive population inversions with the 
third mode, while the first mode holds more than half of the trapped bosons. 



. Figure [3(b) shows the chaotic population dynamics for x 
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Figure 8. [colour] Left (center) figure shows the quantum evolution of the average 
occupations in the three local modes for TV = 30, |Q| = 1, x = 5 (x = —5), /i = 
considering as the initial state the coherent state centered at the initial condition of 
the trajectory in figure [71(a) (figure [JJ(b)). Notice in (b) that the curves for (a\ai) 
and (0303) are almost superimposed. The solid [blue] and dashed [green] curves in (c) 
represent the dynamics of "P su (3) corresponding to the results displayed in (a) and (b) , 
respectively. 



For direct comparison with the previous classical results, the exact quantum 
population dynamics is shown in figure El where we used as initial states the coherent 
states centered at the initial conditions of figure [TJ Figure [HJ (a) displays the quantum 
population dynamics corresponding to the regular trajectory near the SDWS. The 
mean occupation in the second local mode is slightly higher than in the classical 
approximation, but the preferential occupation in the other two modes still prevails. 
Therefore, there is a fair agreement between the exact quantum result and the classical 
approximation when we consider a regular orbit close to the stable SDWS. However, the 
quantum dynamics corresponding to the chaotic trajectory, shown in figure [HI (b), has 
little resemblance to its classical analogue. Unlike the classical approximation, there 
is no preferential occupation of the first mode, which displays values below y. On 
the other hand, the second well does not remain empty, as expected for the irregular 
dynamics near the SDWS. 

The dynamics of < P SU (3) confirms that the classical approximation is more accurate 
for the regular trajectory, as we see in figure [81(c). The state associated with the chaotic 
dynamics has a large purity loss, which corresponds to a large derealization in the phase 
space. Conversely, the time evolution of the state under the regular dynamics preserves 
its similarity to the coherent states, i.e., the state remains well localized in the phase 
space and exhibits relatively high values of purity during its evolution. 
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3.3. Vortex Dynamics 



The vortex states are the only equilibrium points of the classical approximation with 
nonzero imaginary parts. Thus, they are the only ones with nonzero angular momentum 
along the symmetry axis of the trapping potential. To emphasize the imaginary parts of 
the complex variables Wj during their time evolution, we introduce a new set of canonical 
variables: 



yj2N-q\-p{ 



(3.10) 



The real variables qj and pj are directly related to the real and imaginary parts of 
Wj, respectively. Note also that these new dynamical variables must satisfy the condition 
2N > q\ +p1 + q% +p\- We show below only the results for the vortex state parametrized 
by w\ = e 1 ^" = u>2' since the two vortex states differ only by the sense of rotation, as 
previously discussed. 





Figure 9. [colour] Left (right) figure shows the Poincare section at p2 = —y/N/2 for 
N = 30, |fi| = 1, % = — 1 (% = —5) and [i = considering only those trajectories 
with the same energy of the vortex state w\ = e l ~ — u;|. The stable (unstable) 
equilibrium point is represented by the [yellow] circle ([red] cross). 



Figure [9] shows the dynamics in the vicinity of the vortex state for two completely 
different situations. For \ = — 1 and /i = ^ the Poincare section at p2 = —^/N/2 
displays the stable vortex state in a regular region of phase space. However, for \ = —5 
and fj, = yoq, the vortex state is unstable and the dynamics of the system is almost 
completely chaotic. 

According to the equation f)2.15p . the angular momentum of the condensate along 
the symmetry axis of the trapping potential is directly proportional to the su(3) operator 
Js — Ji + <h + J?,- Therefore, the rotational dynamics of the condensate is illustrated 
in Figure [101 where we consider the initial conditions indicated by the [green] star 
markers in figure [91 The solid [blue] (dashed [green]) curve in figure [101 (a) shows the 
classical dynamics of a regular (chaotic) trajectory near the vortex state. Note that 
only the regular trajectory represents a persistent collective rotation of the condensate. 
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Figure 10. [colour] Left (middle) figure exhibits the classical (quantum) dynamics 
of {Js)-, quantity directly proportional to the angular momentum of the condensate 
along the symmetry axis of the trapping potential. The solid [blue] (dashed [green]) 
curve shows the time evolution for N = 30, |0| = 1, x = — 1 (x — ~ 5), \i = and 
the initial condition indicated by the [green] star marker in figure [51(a) (figure 0(b)). 
Right figure shows the dynamics of 'P su (3) for the same initial conditions. 



Therefore, the rotation of the condensate in a preferential sense is present only when 
the vortex state is stable. The chaotic trajectory does not exhibit preferential sense 
of rotation, but presents oscillations of angular momentum with large amplitudes. 
However, the maximum value of (Js(i)} in both trajectories remains bounded by y/SN, 
which corresponds to the mean value of Js for the vortex state, according with (|2.32J) . 

Figure [101 (£») shows the quantum dynamics of rotation, where we consider the 
initial coherent states centered on the initial conditions of figure [10j(a). Note that 
even for N = 30, a relatively small number of trapped bosons, the quantum results 
of (Js(t)) associated with the regular trajectory demonstrate good agreement with the 
classical approximation. However, in the chaotic case, the amplitude of the oscillations 
are strongly attenuated. Thus, as expected, the approximation is not quantitatively 
satisfactory in the chaotic regime, although the absence of a preferential sense of rotation 
is evidenced in both approaches. Figure fTOl(c) displays the dynamics of P su (3) for the 
two orbits previously selected. Note that the purity loss in the chaotic regime is much 
faster and more intense when compared to the regular system, confirming the differences 
found in the precision of the classical approximation. 

4. Conclusion 

In the present work, we have investigated the quantum dynamics of a Bose-Einstein 
condensate confined in a symmetric triple-well potential employing an approximation 
based on the TDVP and the SU(3) coherent states. We have exploited maximally the 
use of coherent states as initial states chosen to be centered in appropriate points of the 
classical manifold, particularly near the equilibrium points of the classical equations of 
motion. This scheme allowed us to capture in the complete quantum calculations some 
of the particularities of the classical solutions. 
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Moreover, in order to evaluate the quality of the classical Hamiltonian evolution, 
as compared to the quantum results, we employed the purity associated with the su(3) 
algebra. The generalized purity measures the departure from the subspace of coherent 
states, in which our classical approximation is restricted, thus enabling us to quantify 
the quality of the approximation. The classical approximation quality depends not only 
on the total number of particles, but also on the location in phase space of the initial 
coherent state and the regularity of the Hamiltonian flow in its vicinity. We concluded 
that even for a relatively small number of particles such as N = 30, which is very 
far from the classical- macroscopic limit, the classical results obtained for the regular 
dynamical regime display a fair agreement with the quantum time evolution, unlike the 
chaotic regime. 

Considering the sub-regime of twin- condensates, we have shown that the self- 
trapping dynamics may be suppressed by the cross-collisions effects on the effective 
tunneling rate, which can become quite important for a large number of trapped bosons. 
Since this sub-regime is classically integrable, its classical approximation exhibits 
excellent results. Whereas, the single depleted well regime, which is characterized by the 
complete depletion of a local mode, can be effectively eliminated in phase space when its 
corresponding fixed points become unstable, giving way to chaos. In the chaotic regime, 
a single Hamiltonian trajectory clearly does not account for the quantum dynamics of 
an initial coherent state at long times. This fact is precisely indicated by the generalized 
purity. In the same line of reasoning, the persistent collective rotation in a preferential 
sense, presented in the vicinity of the vortex states in phase space, can only be properly 
observed when the corresponding fixed points are stable. 
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